Discretized Diffusion Processes. 
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We study the properties of the "Rigid Laplacian" operator, that is we consider solutions of the 
Laplacian equation in the presence of fixed truncation errors. The dynamics of convergence to the 
correct analytical solution displays the presence of a metastable set of numerical solutions, whose 
presence can be related to granularity. We provide some scaling analysis in order to determine the 
value of the exponents characterizing the process. We believe that this prototype model is also 
suitable to provide an explanation of the widespread presence of power-law in social and economic 
system where information and decision diffuse, with errors and delay from agent to agent. 
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Most equations in science are continuous in value. 
While this is a good approximation to reality, natural 
processes often are discrete in microscopic details (e.g. 
atoms and molecules) . We ask the following general ques- 
tions: if the variables are subject to some small threshold 
such that no changes smaller than such a threshold are 
implemented, how do the continuous equations behave? 
Can we properly recover the continuous limit letting the 
threshold vanish? Are there new features emerging in- 
dependent on the threshold value (and therefore robust 
down to the continuous limit)? 

This problem is somewhat general and this is the rea- 
son of its importance. If we consider in particular the 
class of self-organized critical systems then the question 
becomes crucial. We know that for many SOC systems, 
such as the BTW Q| or the Zhang model H, threshold 
and discretisation play a dramatic role in order to avoid 
the diffusive behaviour and to attain the critical one. By 
considering the continuous limit of the equations describ- 
ing the microscopic dynamics of these cellular automata 
and by not preserving a finite threshold we get a trivial 
diffusive behaviour instead of the critical one H . 

In order to investigate this kind of subject, we study 
here one of the oldest and most important equations of 
physics: the Laplace equation V 2 0(x) = 0. Its time 
dependent version, describes the diffusion of heat and 
of particles, and the relaxation of incompressible, non- 
viscous fluids. In imaginary-time it corresponds to the 
Schroedinger equation of quantum mechanics for a free 
particle (the difference vanishes completely if we are in- 
terested in stationary solutions/eigenstates). Therefore, 
we do not overestimate its importance if we consider it to 
be ubiquitous in physics. More recently, it was also con- 
sidered to play a major role in fractal growth processes 
such as Diffusion Limited Aggregation (DLA) [Q] and the 
Dielectric Breakdown Model (DBM) ||. 

It is thus not a surprise that many methods have been 
developed to solve it, at least numerically. In general, 
both time and space are discretized on a lattice. The 
iterative method is the most commonly used method of 



solution, especially for problems such as DLA and DBM 
where the boundary conditions change in time. Essen- 
tially, the method resorts to iterating the equation 
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whose fixed point correspond to the discretized version 
of V 2 4>(x) = 0. Numerically convergence to the fixed 
values <p(i) = (f>(i,t — oo) has to be defined by some 
error: when two successive configurations differ by less 
than some small 5 then the iteration stops. Of course <5 
has to be greater than the machine precision, 5' . This 
means that when the field <fi is of order of unity, it is 
impossible to change its value of quantities less than 5' . 

Therefore it is of extreme importance to explore what 
happens if we explicitcly set a precision in the definition 
of (j), that is, if we say that we can change its value, 
according to (0), only by integer multiples of some chosen 
e. This corresponds to consider the fields 4> 
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Int £ [x] corresponds to taking the integer multiple of e 
closest to x (other possible definitions of the integer part, 
such as the lowest integer closest to x, have been explored 
and have been shown to give the same results). The rele- 
vance of this problem is not only related to the numerical 
solution of Laplacian equation, but it comes also from 
real problems: granular materials, for example, behave 
sometimes as fluids. As such, we can expect them to 
obey some of the laws typical of fluids (hence also to re- 
lax according to the time dependent Laplacian equation) , 
yet their intrinsic granularity forbids movements of quan- 
tities smaller than a single grain. Under this respect we 
refer to the properties of <j) as properties of rigid diffusion 
processes. Discrete-discrete Laplacian fields (f> (discrete 
in space-time, and discrete in the field <p) can also play 
a role in modeling of economic and social systems. In 
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particular, for the simplest case of the Laplacian "har- 
monic" operator the field is computed by averaging over 
the field neighbours. If this is viewed as a social system 
where opinions of players are formed by consulting people 
around, we would like to introduce a more realistic ver- 
sion, where decisions are considered only when a certain 
threshold of information is overcome and the diffusion is 
therefore discretized. This effect known as "rigidity" in 
the literature (^] is indeed responsible of the onset of crit- 
icality ||. Since Self-Organized Criticality is believed to 
play a role both in economic and social system we sug- 
gest that our results can be viewed as a way to describe 
the appearance of criticality in such environments. 

In this Letter we explore the behavior of this Rigid 
Laplacian Model (RLM). As we shall see, non trivial 
properties emerge already in one dimension. 

The standard setup of our simulations is the following: 
we take a Id lattice of L+l sites, and we set the boundary 
conditions 4>(Q) — 1 and cj)(L) — 0. Then we iterate the 
equation 



(i,t + 1) = Int e 
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The iteration stops when there are no more possible rear- 
rangements. The exact solution of V 2 <fi = on the lattice 
with the above boundary conditions is the straight line 
(f>(i) = —if L+l. As we show in Fig.0, the typical solution 
of (H) is instead a quadratic curve. 



FIG. 1. Profile of <f> for a d = 1 system whose length is 
L = 1000. At boundaries 0(0) = 1, = °- 

This can be easily explained by considering that as 
soon as a certain precision s is introduced, the Laplace 
equation can be considered as a first approximation as 



Vd> = ae 
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In the case where the boundary condition are <j)(L) = 
and 4>(0) = 1 we obtain the following equation 
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cf> e (x) = z aeL {j) - ( 1+ 2 aeL ) L +1 ' (5) 



in terms of the reduced variables x' = x/L and e' 



T^aeL 2 one obtains 
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(1 + e'y + 1. 



(6) 



The precise value of a depends on the initial condition 
and it has to be determined with some fitting procedure. 
We can study the dynamics of the system observing that 
the evolution (|l]) always minimizes the functional 
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Eq. (^) maps the RLM into the discrete space-time Gaus- 
sian model. Physically, since 4> represents an electrostatic 
field, E represents the total electrostatic energy of the 
sample, whose minimization in the case of fractal growth 
has been recently studied 0. Letting the field </> take on 
only integer values leads to the discrete Gaussian model. 
The dynamics described by (H|) is therefore aT — dy- 
namics, always decreasing the energy, and the stability 
of solutions different from the constant slope implies the 
presence of local energy minima. Since this process of en- 
ergy minimization stops in local minima, we can perform 
some sort of simulated annealing by perturbing it. We 
increase the value of the field <p in some randomly chosen 
point by ne, then we let the system relax. Indeed, we find 
that the curvature a decreases in time (time is defined in 
terms of perturbation steps; we consider the relaxation 
process to be much faster, as usual in self-organized crit- 
ical models). Each perturbation implies a rearrangement 
of the profile. 

For the energy of the system one has from the eq. (||) 



E cx L 3 a 2 s 2 . 



(8) 



We know that for ordinary diffusion the energy variation 
E is proportional to i~ 3 / 2 . We will show in the conclusion 
that based on that result, in this case one can expect 
E oc t 1 / 2 or more precisely E — Eq oc (t/L 2 ) 1 / 2 , where 
we explicitely introduced a dynamical exponent z for the 
process. Assuming for the moment this result, one has 

EolEq - {t/L z ) 1/2 =>(t/L z ) 1/2 ozE - L 3 a 2 {t/L z ) (9) 

By dcrivating this expression and by requiring that it 
shows no dependence on L (unless for the ratio t/L z ) 
one obtains that the dynamical exponent must be z = 3, 
in agreement with the value z — 3.0 ±0.1 from numerical 
simulations, and 



*(t,L) = a -(t/L z ) 1 / 2 f(t/L z ). 



(10) 



the scaling function f(x) is such that a — > a^, ^ for 
t 3> L z . Such behaviour has been tested in Fig.|| We 
see that for this kind of systems the relaxation process 
to the equilibrium is very slow because of granularity. 
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Indeed, for a continuous Laplacian field the relaxation 
to the stationary state is an exponential process with a 
characteristic time scale to ~ L 2 . 




FIG. 2. Plot in rescaled time of a(t) and a(t) — ao 

Such rearrangements proceed through power-law dis- 
tributed avalanches, whose measure is the number of 
sites that changed their value of <fi. We show in Fig.|| 
the probability p(s) to deal with an avalanche of size s. 
This function can be fitted with a power law of the kind 
p(s) oc s~ T . The value of the exponent r measured from 
the data is 3.8 ± 0.3. To take in account the finite size 
effects of this distribution (the size of an avalanche can- 
not be greater than the size of the system) we can write 
the universal form p(s,L) = s~ T g(s/L), where g(x) is 
a scaling function such that p(s = L) = 0. This scale 
invariant behavior describes the dynamics in a transient 
period where perturbations are accumulated into the sys- 
tem, until a limit profile (not necessarily the straight one) 
is reached. This transient state can be suitably delayed 
by considering a different dynamics where also negative 
perturbations are considered. To describe the evolution 
to the stationary state, we measured the average size (s) 
of the avalanches with respect to the time of the simula- 
tion. 




FIG. 3. Plot of the density function p(s) giving the number 
of times that s sites are updated in an avalanche. 



In Fig]| we present the behaviour of (s(t)) that also 
behaves as a power law with exponent —1/3. A simple 
scaling argument can be presented to explain such be- 
haviour. Indeed, 



(s(t,L)) = Lt- v g(t/L*) 



(11) 



where we explicitely considered the finite size effects in 
the average size. In the above expression, v represents the 
expected exponent and z is the dynamical exponent pre- 
viously introduced and present also in the scaling func- 
tion g. 

The collapse shown in the inset of Fig.|] has been ob- 
tained with the following values: z = 3.0±0.1, v = 
0.34 ± 0.01. The values of these two exponents are re- 
lated, indeed the meaning of the dynamical exponent is 
as usual to determine the scale of time after which the 
process stops. By imposing that at the maximum time 
the average size of the avalanche is independent on L, 
one obtains the scaling relation 
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(12) 



from which v = 1/z = 1/3 as numerically found in the 
limit t « L 3 . 



10 J 



10 




FIG. 4. {s(t)) versus t, in the inset collapse plot for differ- 
ent L's. 

From this scaling relation it's possible to provide an 
argument to determine the value of the exponent r. Since 
the number of avalanches is the same whether they are 
classified according to their sizes (i.e. by means of the 
p(s)ds) or just counted in time (i.e. by means of dt) then 
one obtains p(s) = From the above scaling laws one 
finally obtains 



t = 1 + l/v. 



(13) 



We numerically tested the validity of these relation- 
ships as above mentioned. The "rigidity" of the field 
makes the process sub-diffusive as it is implied by the 
value z = 3 of the dynamical exponent, as opposed to 
the usual z = 2 value. A more intriguing feature of these 
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processes can be obtained by a suitable rescaling of the 
microscopic time scale. As a matter of fact one can check 
the number n(t) of microscopical sweeps needed in order 
to update the values of the field (f> in an avalanche. This 
number n(t) must be proportional to the average size 
of an avalanche at time t and then it varies as t^ 1 ^ 3 , 
as it can be shown by a direct measurement (Fig.|j|). It 
means that the "physical" time of the process defined 
as Y^t'=i n (t') behaves as i 2 / 3 . By rescaling the scale of 
time with this characteristic time scale, the process loses 
any information on the granular nature of the medium 
and the dynamics behaves as a usual diffusive process. 
In this way we demonstrate the initial assumption for 
the energy E in the system. In fact, since it varies for 
ordinary diffusion phenomena as E oc i' -3 / 2 and in this 
case the time of the process t' varies as t' oc t~ x / 3 we 
obtain the above assumption of E oc t 1 !" 1 . 

We also checked the properties of this model in dimen- 
sion larger than d = 1. We found a similar scale invari- 
ant behaviour. In the case of d = 2 we found a value of 
t = 2.65 ± 0.03, for d = 3 we have r = 2.60 ± 0.05 and 
for d = 4 r = 2.55 ± 0.05. 
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FIG. 5. Plot of n(t); the dashed line is y = i 1/3 . 

As described in the introduction, it would be tempting 
then to link this randomness to the known effect of "rigid- 
ity" in Self- Organized systems |^||. Self-Organization 
has often been invoked as a way to describe the invari- 
ance observed in various social and economic systems 
|^| |llj. Indeed we propose the basic mechanism of dif- 
fusion of money and/or information in a random envi- 
ronment (that can stop or reduce the information inter- 
changed between players) as one of the reasons that could 
lead to scale invariance and to avalanche dynamics. As 
regards the economic systems, we believe that this study 
could be used also to describe the time behaviour of a 
quantized field with respect to the value of the trunca- 
tion error. A specific example in this field is be presented 
by considering the spurious effect in the price dynamics 
that could happen by changing the tick at which they are 
traded (i.e. from 1/16 of dollar to 1/100). 



In conclusion, we introduced a model of diffusion, 
where the redistribution process is hindered by some sort 
of rigidity induced by the granularity of the fields. Ran- 
domness is present through a cutoff in the different pos- 
sible values the field can assume. By truncating to a 
certain precision the value of the field a new and rich 
scenario of metastable states appears. The evolution of 
this system through the metastable states is driven by 
an avalanche dynamics with no particular time or length 
scale. By suitably rescaling the microscopic time of the 
evolution, one can use the scaling relations of diffusive 
process to give a theoretical ansatz for the quantities de- 
scribing the system. From computer simulations we find 
a good agreement between numerical data and theoretical 
ones. We believe this sort of process could be responsi- 
ble for the ubiquitous presence of power-law relations in 
the every day life. Clearly our work is readily generalised 
to other systems, notably Naviers- Stokes, KPZ wave 
equations. One may wonder if there is some sort of super- 
universality in their scaling behavior. Indeed a long list 
of similar problems awaits us to study. 
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